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Abstract.  Using  density  functional  theory  with  empirical  van  der  Waals  corrections,  cold  pressure 
curves  were  calculated  and  combined  with  the  quasi-harmonic  approximation  to  study 
thermodynamical  properties  of  several  energetic  molecular  solids.  Vibration  spectra  at  each 
compression  were  calculated  and  used  for  including  temperature  and  zero-point  energy  contributions  to 
the  free  energy.  Equilibrium  properties  at  temperatures  of  experiments,  as  well  as  hydrostatic  equations 
of  state,  specific  heat  capacities,  and  coefficients  of  thermal  expansion,  were  obtained  and  compared  to 
experiment. 
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INTRODUCTION 

To  perform  meso-scale  and  continuum-level 
simulations  of  detonating  energetic  materials 
(EMs),  accurate  mechanical  and  thermal  properties 
are  required  to  provide  meaningful  results  on 
explosive  performance,  as  well  as  materials 
response  to  a  variety  of  conditions  [1]. 
Unfortunately,  experiment  falls  short  of  providing 
these  properties  over  a  wide  range  of  pressures  and 
temperatures,  and  is  hindered  by  inaccuracies  from 
shear  introduced  by  the  compressive  media  used  in 
diamond  anvil  cells  [2].  As  such,  theoretical 
avenues  for  the  acquisition  of  key  input  parameters 
for  large-scale  modeling  of  EMs  have  been  a  major 
thrust  for  the  EM  modeling  community. 

Initial  efforts  focused  on  calculating  equations 
of  state  (EOSs)  and  thermal  properties  using  pure 
density  functional  theory  (DFT)  methods;  however, 
DFT  with  generalized  gradient  approximation 
(GGA)  greatly  over-predicted  equilibrium  volumes 


by  an  average  of  10%  compared  to  experiment  due 
to  inadequate  description  of  dispersive  van  der 
Waals  (vdW)  interactions  [3,4].  More  recently, 
attempts  to  economically  incorporate  vdW  into 
DFT  (DFT+vdW)  led  to  various  treatments  using 
pair-wise  C6r"6  potentials  that  are  simply  added  to 
DFT  at  large  distances  [5,6].  While  these 
approaches  have  been  shown  to  vastly  improve  the 
prediction  of  equilibrium  volumes  in  EMs,  they 
systematically  under-predict  these  volumes  by 
about  3%  [7].  Fortunately,  this  under-prediction 
was  shown  to  be  due  to  the  exclusion  of  thermal 
and  zero-point  energy  (ZPE)  effects  on  the 
crystalline  environment  [8].  By  including  vdW, 
thermal,  and  ZPE  effects  into  DFT  (DFT+vdW+T), 
equilibrium  volumes  were  predicted  within  1%  of 
those  obtained  by  experiment,  and  very  accurate 
EOSs  were  obtained  [8]. 

In  addition  to  obtaining  mechanical  properties 
of  EMs  by  DFT+vdW+T,  thermal  properties  are 
also  of  substantial  interest.  Therefore,  the  goal  of 
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this  work  is  to  calculate  thermal  properties,  such  as 
heat  capacities  and  coefficients  of  thermal 
expansion  (CTEs),  for  EMs  using  the 
DFT+vdW+T  approach,  and  compare  these 
predicted  properties  to  those  of  experiments. 

COMPUTATIONAL  DETAILS 


The  DFT+vdW+T  approach  used  in  this  work  is 
based  on  standard  DFT  with  the  Perdew-Burke- 
Ernzerhof  exchange-correlation  functional  and  the 
empirical  vdW  treatment  of  Neumann  and  Perrin 
[6].  The  vdW  method  employed  uses  C6f6  pair¬ 
wise  potentials  multiplied  by  a  damping  function 
that  prevents  divergence  of  1/r6  at  small  r. 


dJ>')= 


1  -  exp 


D 

\^ab  y 


(1) 


v  l  jy 

The  damping  function,  dAB  in  (1),  contains  two 
parameters  that  were  fit  to  experimental  data  for 
cold  (T  <  5 OK)  molecular  crystals:  the  form  factor 
n,  which  alters  the  shape  of  the  damping  function; 
and  the  crossover  distance  RAB,  which  defines  the 
short-distance  cutoff  for  vdW  interactions. 

The  DFT+vdW  approach  was  used  to  calculate 
hydrostatically  compressed  and  expanded  cold- 
pressure  structures  for  which  the  atoms  were 
relaxed  to  a  maximum  force  of  0.03  eV/A.  For 
each  compression/expansion,  the  dynamical  matrix 
is  calculated  by  finite  displacement  with  a  step  size 
of  ±0.05  A.  Diagonalizing  these  matrices  give 
vibration  spectra  under  the  quasi-harmonic 
approximation,  which  can  then  be  used  to  calculate 
free  energies  and  thus  equilibrium  volumes  for  the 
EMs  studied  at  various  temperatures;  see  Ref.  [8]. 
Additionally,  the  vibration  spectra  were  used  to 
calculate  heat  capacities  and  CTEs  up  to  melt 
temperatures  for  comparison  to  experiment. 


RESULTS  AND  DISCUSSION 


PETN-I  determined  by  experiment  [9]  at  300  K  is 
589.5  A3,  and  is  predicted  to  be  590.7  A3  by 
DFT±vdW±T.  Likewise,  NM  at  4.2  K  was 
determined  [10]  to  be  275.3  A3,  but  was  also 
predicted  to  be  275.3  A3  by  DFT±vdW±T. 
Additionally,  the  resulting  bulk  moduli  for  NM  and 
PETN-I  are  7.9  and  12.8  GPa  respectively,  which 
compares  favorably  to  the  values  8.3  and  12.3  GPa 
obtained  by  experiment  [11,12].  This  demonstrates 
the  ability  of  DFT±vdW±T  to  reasonably  predict 
bulk  moduli. 
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Figure  1.  EOSs  for  NM  and  PETN-I  showing 
agreement  between  experimental  data  (triangles)  and 
data  calculated  using  DFT±vdW±T  (crosses).  Solid  line 
is  B-M  fit  for  calculated  data. 
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As  seen  in  Figure  1,  EOSs  were  calculated  for 
both  nitromethane  (NM)  and  PETN-I  using  pure 
DFT,  DFT±vdW,  and  DFT±vdW±T.  Evidenced  by 
their  relation  to  experiment  (triangles)  pure  DFT 
over-predicts,  while  DFT±vdW  under-predicts  the 
EOSs.  Only  when  temperature  and  ZPE  effects  are 
included  is  good  agreement  with  experiment 
obtained.  For  instance,  the  unit  cell  volume  for 


The  predictive  value  of  the  DFT±vdW±T 
approach  for  mechanical  properties  of  EMs 
suggests  its  potential  for  predicting  thermal 
properties  such  as  heat  capacities  or  CTEs  over  a 
wide  range  of  pressures  and  temperatures.  In  fact, 
as  seen  in  Figure  2  (top  and  middle),  fairly  good 
agreement  can  be  seen  in  the  cases  of  NM  and 
PETN-I  at  ambient  pressure;  however,  the  values 
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obtained  from  calculation  diverge  from  those 
obtained  from  experiment  approaching  the  melt 
temperature.  Additionally,  the  bottom  plot  in 
Figure  2  shows  the  rather  poor  agreement  of  the 
calculated  CTE  as  compared  to  that  obtained  by 
experiment  for  PETN;  where  other  EMs  showed 
even  worse  agreement,  and  thus  are  not  reported 
here  for  compactness.  Clearly  the  small 
discrepancies  in  heat  capacities,  and  the  large 
discrepancies  in  CTEs,  result  from  inaccuracies  in 
the  vibration  spectrum,  particularly  the  low 
frequency  modes,  to  which  heat  capacity  and  CTE 
are  most  sensitive,  and  free  energy  is  least 
sensitive.  To  sufficiently  sample  intermolecular 
modes,  a  1x1x2  cell  was  used  for  PETN-I  such  that 
the  minimum  dimension  of  the  most  compressed 
structure  was  not  less  than  8.5  A.  Likewise,  a 
2x2x1  cell  was  used  for  NM. 

To  rule  out  both  inaccuracies  in  forces,  and 
sampling  of  the  dynamical  matrix  beyond  the 
harmonic  limit,  the  force  convergence  criteria  was 
lowered  to  0.01  eV/A,  while  the  finite 
displacement  step  size  was  lowered  to  ±0.02  A. 
These  changes  resulted  in  slight  improvement  to 
heat  capacity,  and  significant  improvement  to  CTE 
for  PETN-I  as  seen  in  Figure  2  (middle  and 
bottom).  Nonetheless,  the  goal  of  achieving 
excellent  quantitative  agreement  with  experiment 
had  not  been  achieved.  Additionally,  it  was  found 
that  despite  the  lack  of  a  polymorphic  change,  the 
calculated  low-frequency,  intermolecular  modes 
exhibited  non-monotonic  behavior  with  volume, 
suggesting  that  the  problem  may  be  due  to  errors 
introduced  by  the  vdW  damping  function  at 
intermolecular  distances.  As  a  final  note,  the  error 
for  the  calculated  equilibrium  volume  of  PETN-I 
increased  from  ±0.2%  to  ±1.2%,  which  can  be 
explained  by  the  fact  that  the  damping  function 
was  fit  to  lattice  parameters  of  low-temperature 
crystals  with  ZPE  effects  still  included. 

To  get  a  sense  of  how  the  vdW  damping 
function  might  affect  the  calculation  of  the 
dynamical  matrix,  the  vdW  interaction  for  two  N 
atoms  was  plotted  along  with  its  first  and  second 
derivatives  with  respect  to  distance;  see  Figure  3. 
The  narrow  peak  in  the  first  derivative,  and 
subsequently  the  resulting  rapid  change  in  the 
second  derivative,  can  greatly  affect  interactions 
within  the  dynamical  matrix  for  atoms  separated  by 
distances  close  to  the  crossover  distance  when  both 


positive  and  negative  displacements  are  made  for 
the  same  atom.  Further,  the  bulk  of  intermolecular 
interactions  for  EMs  occur  within  ±2  A  of  the 
crossover  distance,  undoubtedly  leading  to  poor 
description  of  intermolecular  modes  with  respect  to 
volume. 


Figure  2.  Specific  heat  capacity  for  NM  (top)  and 
PETN-I  (middle),  as  well  as  CTE  for  PETN-I  (bottom). 
For  PETN-I,  the  high  accuracy  curves  correspond  to 
lowered  force  convergence  criteria  and  smaller  step  size 
for  finite  displacement. 

CONCLUSIONS 

In  summary,  DFT±vdW±T  predicts  well  the 
EOSs  of  EMs  over  ranges  available  from 
experiment,  but  does  only  moderately  well  for  heat 
capacities  and  poorly  for  CTEs  -  both  being 
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sensitive  to  low  frequency  modes.  The 
discrepancies  seem  to  be  due  to  inaccuracies  in 
low-frequency,  intermolecular  modes  brought  on 
by  the  non-monotonic  behavior  of  the  second 
derivative  of  the  empirical  vdW  interaction.  This  is 
a  direct  effect  of  the  damping  function  form  used  to 
shut-off  vdW  interactions  at  small  distances. 
Further  refinement  to  the  damping  function  form, 
as  well  as  fits  to  crystallographic  data  extrapolated 
to  cold  pressure  values  to  redundancy  in  the 
treatment  of  ZPE  effects,  should  greatly  improve 
accuracy  in  the  prediction  of  thermal  properties  of 
EMs,  and  further  refine  accuracy  for  EOSs. 


0.0  1.0  2.0  3.0  4.0  5.0  6.0  o  7.0  8.0 

N-N  Interatomic  Distance  (A) 

Figure  3.  Calculated  vdW  binding  energy  curve  (top) 
for  prototypical  two  nitrogen  atom  system  with  first 
spatial  derivative  (middle)  and  second  spatial  derivative 
(bottom)  demonstrating  inaccuracies  introduced  near 
the  crossover  distance  (Rn-n  =  3.384  A). 
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